#!/usr/bin/perl
# 2014/01/07 From capillary sequence results, split sequence for alignment

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use File::Basename;
use Carp qw(croak);
$| = 0;

if($ARGV[0] eq '--help'){
	die "Usage: perl $0 -d [directory or data]\n";
}



my %opt;
getopts("d:", \%opt);
my $DIR = $opt{d};

my $FILE_FASTA = $DIR . '/' . 'original.fasta';
my $FILE_blast = $DIR . '/' . 'blast_result.txt';
my $FILE_fastq1 = $DIR . '/' . 'left.fasta';
my $FILE_fastq2 = $DIR . '/' . 'right.fasta';
my $FILE_table = $DIR . '/' .  'table.txt';
my $LINKER_LENGTH = 34;	# Both ChIA-PET and HiD are same

#---------------------------------------
# making FASTA file
#---------------------------------------
my $fh_out_fasta = IO::File->new($FILE_FASTA, 'w') or die "cannot write $FILE_FASTA: $!";
my $fh_out_table = IO::File->new($FILE_table, 'w') or die "cannot write $FILE_table: $!";
opendir DIR, $DIR or die "cannot open file: $!";
while(my $file = readdir DIR){
	if($file =~ /(.+)\.seq$/){
		my $ID = $1;
		$fh_out_fasta->print(">$1\n");
		$fh_out_table->print("$ID\t");
		my $fh_in = IO::File->new($file) or die "cannot open $file: $!";
		while($_ = $fh_in->getline()){
			s/\r?\n//;
			$fh_out_fasta->print("$_");
			$fh_out_table->print("$_");
		}
		$fh_out_fasta->print("\n");
		$fh_out_table->print("\n");
		$fh_in->close();
	}
}
$fh_out_fasta->close();
$fh_out_table->close();

#---------------------------------------
# BLAST
#---------------------------------------
my $blast_CMD =<<"EOF";
blastn -query $FILE_FASTA -strand 'plus' -db HiD_linker -task blastn -out $FILE_blast -outfmt "6 qseqid sseqid qstart qend sstart send" -max_target_seqs 1;
EOF

### example output
#ID_0	ChIA-PET_linkerBB	24	57	1	34
#ID_1	ChIA-PET_linkerBB	24	57	1	34
#ID_2	ChIA-PET_linkerBB	24	57	1	34
#ID_3	ChIA-PET_linkerBB	26	57	3	34
#ID_4	ChIA-PET_linkerBB	24	57	1	34

system($blast_CMD) == 0 or die "BLAST failed; $!";




#---------------------------------------
# split file
#---------------------------------------
{
	my $fh_blast = IO::File->new($FILE_blast) or die "cannot open $FILE_blast: $!";
	my $fh_table = IO::File->new($FILE_table) or die "cannot open $FILE_table: $!";
	my $fh_seq1 = IO::File->new($FILE_fastq1, 'w') or die "cannot write $FILE_fastq1: $!";
	my $fh_seq2 = IO::File->new($FILE_fastq2, 'w') or die "cannot write $FILE_fastq2: $!";

	while(my $blast = $fh_blast->getline()){
		my $table = $fh_table->getline();
		$blast =~ s/\r?\n//;
		$table =~ s/\r?\n//;

		my ($bid, $linker, $qstart, $qend, $sstart, $send) = split /\t/, $blast;
		my ($tid, $seq) = split /\t/, $table;

		### sometimes blast cannot find linker sequnece and skipped those read
		### then skip treatment of those read and move to next
		while($tid ne $bid){
			$table = $fh_table->getline();
			($tid, $seq) = split /\t/, $table;
			if($tid eq ''){
				last;
			}
		}

		if($bid ne $tid){
			die ("blast and table id not matched ! at \n$blast\n$table\n");
		}

		# check barcode sequence
		my $barcode = 'NA';
		my $barcode_target_seq = substr($seq, $qend);
		if($barcode_target_seq =~ m/GAC([ATGC]{5})GCA/){
			$barcode = $1;
		}

		# combin id
		my $ID = '@' . $bid . ':' . $barcode;

		# adjusting the alingment position in case of partial alignment
		if($sstart != 1){
			$qstart = $qstart - ($sstart - 1);
		}
		if($send != $LINKER_LENGTH){
			$qend = $qend + $LINKER_LENGTH - $send;
		}

		# start begins from 0
		$qstart -= 1;
		$qend -= 1;

		### frist half part
		my ($seq_1, $seq_2);
		my ($quality_1, $quality_2);
		{
			my $start = $qstart - 25;
			if($start < 0){
				$start = 0;
			}
			$seq_1 = substr($seq, $start, 25);
			$fh_seq1->printf("%s\n%s\n", $ID, $seq_1);
		}
		{
			my $start = $qend + 1;
			$seq_2 = substr($seq, $start, 25);
			$fh_seq2->printf("%s\n%s\n", $ID, $seq_2);
		}

	}
	$fh_blast->close();
	$fh_table->close();
	$fh_seq1->close();
	$fh_seq2->close();
}


